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Abstract 

A method is proposed for generating compact fractal disordered media, by generalizing the ran- 
dom midpoint displacement algorithm. The obtained structures are invasive stochastic fractals, 
with the Hurst exponent varying as a continuous parameter, as opposed to lacunar deterministic 
fractals, such as the Menger sponge. By employing the Detrending Moving Average algorithm 
[Phys. Rev. E 76, 056703 (2007)], the Hurst exponent of the generated structure can be subse- 
quently checked. The fractality of such a structure is referred to a property defined over a three 
dimensional topology rather than to the topology itself. Consequently, in this framework, the Hurst 
exponent should be intended as an estimator of compactness rather than of roughness. Applica- 
tions can be envisaged for simulating and quantifying complex systems characterized by self-similar 
heterogeneity across space. For example, exploitation areas range from the design and control of 
multifunctional self-assembled artificial nano and micro structures, to the analysis and modelling 
of complex pattern formation in biology, environmental sciences, geomorphological sciences, etc. 
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I. INTRODUCTION 



Real world materials, like stable and metastable liquids, glasses, defective crystals, and 
structures that evolve from non equilibrium processes, always exhibit a certain degree of dis- 
order. On the other hand, the design of artificial heterogeneous structures and the accurate 
control of their structural, chemical or orientational disorder is an area of intensive investiga- 
tion for fundamental and technological interest [Ill21l3lSll51El[ZllSllSlliniinilI2llI3]. Recent 
developments in nanophotonics have shown that it is possible to make use of the intrinsic 
disorder in natural and photonic materials to create new optical structures and function- 
alities . Fractal bodies, such as Menger sponges [15], have demonstrated the ability to 
perform novel functions as localize electromagnetic and acoustic waves [HI H7J HE] or en- 
hance super-liquid repellency or dewettability [191 [20]. The reconstruction of heterogeneous 
media, from the knowledge of the correlation function, is a challenging inverse problem. 
Any reconstruction procedure requires to be effectively controlled in order to enable the 
prediction, design and implementation of structures exhibiting the desired electromagnetic, 
transport or biological functions [2HE21I2S1I21ES1I2S1I2ZII2B]. 

The disorder degree of a medium can be quantified in terms of the two-point correlation 
function C(ri,r 2 ) =< /(ri)/(r 2 ) > of a relevant quantity /(r) (e.g. dielectric function, 
porosity, density), where ri,r-2 are two arbitrary points in the system. For statistically 
isotropic media, C(r 1 ,r 2 ) depends only on the distance A = ||r 2 — rx|| between two points, 
thus is written as C(A). For fully uncorrelated systems, like the ideal gas, the correlation 
function is a simple exponential, C(A) oc exp(— A/a), while for fully ordered media, like 
the perfect lattice, the correlation function C(A) is a constant. Generally, heterogeneous 
materials exhibit forms of correlation which are intermediate between those of the ideal gas 
and the perfect lattice. 

Fractional Brownian functions are characterized by a correlation function depending as 
a power-law on A [2HJ ED] • Such a correlation may reasonably account for the intermediate 
behavior, between the fully uncorrelated exponential and the fully correlated constant de- 
cay, exhibited by real disordered media. The power-law correlation of fractional Brownian 
functions can be expressed by the power-law decay of the variance: 

<[Mr + A)-Mr)] 2 ) =al\\\\\ 2H , (1) 

with f H (r) : R d -> R, with r = (x ± , x 2 , x d ), A = (A x , A 2 , X d ) and ||A|| = 

2 



yAi + X\ + ... + . H is the Hurst exponent, which is related to the fractal dimension 
through the relation D = d + 1 — H, with d the embedded Euclidean dimension. H ranges 
from to 1, taking the values H = 0.5, H > 0.5 and H < 0.5 respectively for uncorre- 
cted, correlated and anticorrelated Brownian functions. Concepts as scaling, criticality and 
fractality have been proven useful to model dynamic processes [3Tj : stress induced morpho- 
logical transformation [32J; isotropic and anisotropic fracture surfaces [331 Ell ESI [361 137] : 
static friction between materials dominated by hard core interactions |38j; elastic and contact 
properties [321 SHI SH H2J SSI SI] ; diffusion and transport in porous and composite materials 
[121 SSI S3 SB]; mass fractal features in wet/dried gels, liposome and colloids [1HJ EH 1ST] ; 
physiological organs (e.g. lung) [52]; polarizabilities [53], hydrophobicity of surfaces with 
hierarchic structure undergoing natural selection mechanism [54J and solubility of nanopar- 
ticles [55]. Several quantification methods have been proposed to accomplish accurate and 
fast estimates of power-law correlations at different scales [561 EU EHl EH [601 ED E2l [631 El] ■ 

In the present work: 

1. A compact invasive stochastic fractal structure is obtained by generalizing the random 
midpoint displacement (RMD) algorithm to high-dimension (Section II). In topological 
dimension d — 3, a fractal cube with size N\ x N2 x A^3 is obtained, whose relevant 
property (e.g. density, dielectric function, porosity) is described as a three-dimensional 
fractional Brownian field As opposed to deterministic lacunar fractals, such as 
the Menger sponge whose fractal dimension is equal to 2.7268, the obtained medium is 
a stochastic invasive fractal. By varying the Hurst exponent, which is given as input of 
the algorithm, between and 1, the fractal dimension D is continuously varied between 
D = 3 and D = 4. It is worthy of remark that in this case the fractal dimension D 
refers to the compactness rather than to roughness as in the case of fractal surfaces 
and interfaces. Such a structure can be used for describing complex materials with 
correlation function exhibiting a power-law dependence over distance and arbitrary 
degree of heterogeneity expressed by H. 

2. The degree of disorder of the medium is quantified in terms of the Hurst exponent, 
whose estimate is provided by the Detrended Moving Average (DMA) algorithm [HS] 
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(Section III). The algorithm calculates the generalized variance <tdma(s) of the frac- 
tional Brownian fields /ff(r) around the three-dimensional moving average function 
fuij'^s). The generalized variance cfdma{s) is estimated over sub-cubes with size 
«i x «2 x ris and, then, summed over the whole domain Ni x N 2 x N3 of the fractal 
cube. The value of udma(s) for each sub-cube is then plotted in log-log scale as a 
function of s = \/nf + n 2 + nj$. The linearity of this plot guarantees the power-law 
dependence of the correlation over the investigated range of scales. The slope of the 
plot yields the Hurst exponent H of the fractional Brownian field. 

II. GENERATION OF COMPLEX HETEROGENEOUS MEDIA 

The random midpoint displacement algorithm is a recursive technique widely used for 
generating fractal series and surfaces. In d — 1, a fractional Brownian walk is obtained 
starting from a line of length N. At each iteration j, the value at the midpoint is calculated 
as the average of the two endpoints plus a correction which scales as the inverse of the length 
of half segment with exponent H. The fractal dimension is D = 2 — H, varying between 
1 and 2. Fractal surfaces with desired roughness can be generated starting from a plane, 
whose domain is a regularly spaced square lattice with size N 1 x N 2 , %\ = 1,2,..., N 1 and 
i 2 = 1,2,..., N 2 . First, the square is divided in four subsquares. Then, the value in the center 
of the square is calculated as the average of the values at the four vertices plus a random 
term. Then, the four values at the midpoint of the edges are obtained as the average of the 
values at the two adjacent vertices plus a random term. The process is repeated until the 
fractal surface is obtained. The fractal dimension is D = 3 — H, varying between 2 and 3. 

Here, a high-dimensional implementation of the algorithm is presented for obtaining a 
compact body, exhibiting power-law correlation over a wide span of scales. 

A. d- dimensional random midpoint displacement algorithm 

Here, the random midpoint displacement algorithm is generalized to be operated on arrays 
with arbitrary euclidean dimension d. The procedure is generalized by defining the function: 




(2) 



k 
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which is defined at the center of a hypercubic equally spaced lattice. The sum is calculated 
over the k endpoints of the lattice. The quantity <jj^ is a random variable defined at each 
iteration j whose explicit expression is worked out below. We start from Eq. ([I]), that is 
written as: ^ 

([Mr + A) - f H {r)f) = ol + A3 + ... + Ajj) , (3) 

by considering a hypercubic lattice with size Ni = N 2 = ... = N d = N, Eq. (J3]) writes as: 

<[/^(r + iV)-^(r)] 2 > = a 2 (7div) . (4) 
By dividing each lattice size by a factor 2 at each iteration, Eq. Q is rewritten: 



2.H 

[Mr + 2j) - Mr)?} =^(^1 • (5) 

Moreover, at each iteration, the following relation holds: 

[f H (r + |) - f H (r)]^ = i (W + ^) - Mr)] 2 ^ + a 2 , . (6) 

because the value of /ff(r) at each step of the RMD algorithm is calculated over two points 
for d = 1 (extremes of a segment), four points for d = 2 (vertices of a square), eight points 
for d = 3 and so on. By using Eq. ([5]), Eq. ^ becomes: 

2H „ / ^ \ 2H 



VdN\ _ a 2 a ( VdN\ 



^ I ] = H I ^rr ) + > ( 7 ) 



and the term cr| d writes: 



,_ x 2H 

VdN 



°l = ° 2 o ^ [1-2^)] , (8) 



2^ 

which is the generalized form of the relationship holding for d = 1 [65] . 



B. Three-dimensional media 



Here, Eqs. ( 2p ) are used to yield a compact random fractal with embedded Euclidean 



dimension d = 3. As already stated, the relevant property of the disordered medium is 
defined by the scalar function /jf(r) : IR 3 — > IR, with r = (xi, X2, X3). To generate the 
fractal structure, a cubic array with size Ni x N 2 x is defined. Initially, the cube is fully 
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FIG. 1: (Color online) Scheme of the elementary steps of the implementation described in Section 
II. Step 1: the homogeneous cube is divided in 27 subcubes (a). Step 2: the function is calculated 
at the subcubes located at the eight vertices (b). Step 3: the function is calculated at the subcube 
located in the center (c). Step 4: the function is calculated at subcubes located at the center of 
the six faces (d). Step 5: the function is calculated at the subcubes located at the midpoint of the 
12 edges. 

homogeneous with the function describing the fractal property of the medium is taken as a 
constant, e.g. fuij") = 0. 

Then, the algorithm is implemented at each iteration j according to the following steps: 

Step 1. The cube is divided in 27 subcubes, as shown in Fig. [IJa). 

Step 2. The values of the function /#j(r) are seeded as random variables of average value 
and standard deviation <jj 3 at the eight subcubes located at the eight vertices, as 
shown in Fig. [l^b). 

Step 3. The values of the function /^ (r), at the central sub-cube, are calculated as follows: 

1 8 

k=l 

where the sum is performed over the /c-values taken by the function f(r) at the sub- 
cubes located at the eight vertices of the main cube. The random variable 0^3 is 
denned as: ^ 

*h = °l(^J [1"2 2 ^] ■ (10) 

Step 4. The values of the function //rj(r), defined over the subcubes located at the center 
of the six faces, are calculated as follows: 

1 4 

^•( r ) = iE^ r )+^ > en) 
fc=i 



where the sum is performed over the k- values taken by the function fHj(f) & t the 
subcubes placed at the four vertices of each face. The random variable 0^2 is defined 
as: 



2 2 

a 3,2 = °o 



V2N 



2H 



[1 - 2 2 ( H -V] 



(12) 



Step 5. The values of the function f'H,j{ r ) at the subcubes located at the midpoint of the 
twelve edges, are calculated according to: 

1 2 

/Hj(0 = 2£/*( r ) + *i,i , (13) 

k=l 

where the sum is performed over the fc-values taken by the function f(r) at the sub- 
cubes located at the endpoints of the 12 edges. The random variable 0^1 is defined 
as: 

Hence, the function ///j (r) at the subcubes located at the midpoint of each edge takes 
a value given by the average of the function at the endpoint subcubes plus the random 
variable a^\. 





(a) (b) (c) 

FIG. 2: (Color online) Fractal media generated according to the procedure reported in Section II. 
The Hurst exponents are respectively H = 0.2 (a), H = 0.5 (b), H = 0.8 (c). The color map is 
rescaled in order that white dots correspond to the minimum values and darkest blue dots to the 
maximum values of /j?(r). 



The first run of the routine results in 27 subcubes, characterized by the values of the 
function /#(r) described above. The structure obtained at the first iteration is shown in 
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(a) (b) (c) 

FIG. 3: (Color online) Granular structures obtained by considering spheres with the radii described 
by the function /j/(r). The Hurst exponents are respectively H = 0.2 (a), H = 0.5 (b) and H = 0.8 
(c). 

Fig. [l](e) . The steps 1-5 are iteratively repeated for each of the 27 subcubes. At the second 
run, a number of subcubes equal to 27 x 27 is obtained. Eventually, the number of subcubes 
will be equal to (3 J ") d , where j is the iteration number and d = 3. 

In Fig. D the fractal cubes with H = 0.2 (a), H = 0.5 (b) and H = 0.8 (c) are shown 
at iteration j = 9. The Hurst exponent H is the input of the generator, which determines 
the heterogeneity of the resulting microstructure. The color map is rescaled in such a way 
that the lightest cubes corresponds to the value of at the initial stage of the medium, 

darker dots to the maximum values obtained by applying the routine to the function /fr(r). 

An alternative representation can be obtained by relating the fractional Brownian func- 
tion to the grain size. Granular structures with different Hurst exponents are shown 
in Fig. [3j In this case, since the fractal property is referred to the size of the grain, whose 
shape is spherical. One can notice that small and large grains are located together for an- 
ticorrelated cubes with H = 0.2. Conversely, as the Hurst exponent increases, the grains 
segregate according to the size. Smaller grains are more likely to be grouped with small 
grains and viceversa. In Fig. |4j granular structures with the same Hurst exponent (H=0.5) 
and different average grain size are also shown. 
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(a) (b) (c) 

FIG. 4: (Color online) Granular structures obtained by considering spheres with the radii described 
by the function /#(r) with three different average grain size. The Hurst exponent is H = 0.5 for 
all the three cases. 

III. ESTIMATING THE HURST EXPONENT OF RANDOM FRACTAL MEDIA 

The aim of this section is to implement an independent measure of the Hurst exponent 
of the heterogeneous structure generated according to the procedure in Section II (shown 
in Figs.[2j[3j[4j The detrending moving average algorithm has been proposed in [56] for 
estimating the Hurst exponent of high-dimensional (d > 2) fractals. In the present work, 
the algorithm will be implemented on the heterogeneous structures generated in Section II. 

The core of the algorithm is the generalized variance crjj MA (s), that for a three- 
dimensional structure, is written: 



4ma( S ) = y Yl [f H ^ ~ fni,n 2 ,n 3 (r) , (15) 
V 



where /ij(r) = /ff(xi, #2; ^3) is the fractional Brownian field with i\ = 1,2,..., N, 12 
1,2, N, i 3 = 1,2, N. The function f nil n 2 ,n 3 (r) is given by: 



fm. 



n 2 ,n 3 \ 



£2, x 3 



(16) 



fel k 2 k :i 

with the size of the subcubes (n 1; n 2 , n 3 ) ranging from (3,3,3) to the maximum values 
( n imax, n2max, n 3max) ■ v = n i n 2n3 is the volume of the subcubes. The quantity V = 
(JVi — ni max )(N2 — n,2max)(N 3 — nzmax) is the volume of the fractal cube over which the 



averages / are defined. As observed in 156] , Eqs. ( 15 ) and ( 16 ) are defined for any geometry of 



the sub-arrays. However, in practice, sub-cubes with n\ = ri2 = n% are computationally more 
suitable for avoiding spurious directionality and biases in the calculations. The generalized 
variance o~ 2 DMA (s) scales as {n\+n\+n\ ) H because of the fundamental property of fractional 
Brownian functions Eq. 0. 



Eqs. (15) and (16) correspond to the isotropic implementation of the algorithm. The 
isotropy follows from the definition of the average /m^^ W, which is obtained by summing 
all the values taken by at the subcubes centered in r. As explained in [56], the 

implementation can be made anisotropic for fractals having a preferential growth direction, 
as for example biological tissues, epitaxial layers, crack propagation. The anisotropy is 
accomplished by varying the sum indexes according to mi = int(ni0i),m2 = int (7^2^2)5 
rris = int(n3#3). Upon variation of the parameters 81, 82, 83 in the range [0, 1] the indexes 



of the sums in Eqs. (15) and (16) are set within each subcube. In particular, r coincides 
respectively with (a) one of the vertices for 81, 82, 83 = or 1 or with (b) the center for 8\ = 
82 — 83 = 1/2. The values 81 = 8 2 = 83 = 1/2 correspond to the isotropic implementation, 
while 8\ = 8 2 = 83 = and 8\ = 8 2 = 83 = 1 correspond to the directed implementation. 



In d = 3, the isotropic implementation implies that the function defined by Eq. (16) is 
calculated over subcubes whose center is r. Conversely, the anisotropic implementation 
implies that the function f ni ,n 2 ,n 3 (r) is calculated over subcubes having one of the eight 
vertices coinciding with r . 

In order to calculate the Hurst exponent of the heterogeneous structure, the algorithm 
is implemented through the following steps. The function f ni ,n 2 ,n 3 ( r ) is calculated over 
different subcubes, by varying n\, n 2 , n 3 from 3 x 3 x 3 to ri\ max x ri2max x n 3max . The 
maximum values n\ max , ri2max , n 3max depend on the size of the whole fractal. To minimize 
finite-size effects, it should be ri\ max « Ni, n 2ma:r « N 2 , n 3max « N 3 . The next step 



is the calculation of the difference fuij') — /ni,n 2 ,n 3 ( r ) in parentheses of Eq. (15) for each 
sub-cube ri\ x n 2 x n 3 . For each subcube, the corresponding value of o~ 2 DMA (s) is calculated 
and finally plotted on log- log axes. The log-log plot of o~ 2 DMA (s) as a function of s , yields a 
straight line with slope H, on account of the following relationship: 

° 2 dma(s) ~ K + n l + n D H ~ s " ■ ( 17 ) 

In Fig. [5J the log- log plots of o~ 2 DMA (s) vs. s are shown for fractal cubes generated according 
to the procedure described in Section II. The cubes have Hurst exponents ranging from 0.1 
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FIG. 5: (Color online) Log-log plot of & DMA (s) for fractal media with size Ni x N2 X N3 = 
1025 x 1025 x 1025. The fractal media are generated by the algorithm proposed in Section II, with 
Hurst exponent varying from 0.1 to 0.9 with step 0.1. Dashed lines represent the linear fits. The 
Hurst exponents can be estimated by means of the 3d-DMA algorithm as explained in Section III. 

to 0.9 with step 0.1 and size 1025 x 1025 x 1025. Dashed lines represent the linear fits, 
with errors and linear regression coefficients shown in Table I. The plots of & DMA (s) as a 



function s are linear according to the power-law behavior expected on the basis of Eq. (17). 



Deviations from the full linearity can be observed particularly at the extremes of the scale. 



IV. CONCLUSIONS 



We have put forward an algorithm to generate a fully compact heterogeneous medium, 
whose fractal dimension can be continuously varied upon varying the Hurst exponent be- 
tween and 1. The generation method is based on a generalization of the random midpoint 
displacement algorithm. In order to check the accuracy of the generator, the Hurst ex- 



ponent of the fractals can be estimated by using the high-dimensional variance a 



DMA 



» 



defined by Eq. (15). We envisage fruitful applications, e.g. in three-dimensional medical 



image analysis, where density or granularity patterns could be interpreted synthetically by 



11 



TABLE I: Hurst exponents Hi, H2, H3, and linear regression coefficients pi, p2, P3, Pa of 
curves like those of Fig. [5] obtained as average of 10 realizations of heterogeneous media with 
Ni x N2 x A3 = 1025 x 1025 x 1025. Hi and H2 have been calculated by linear fit over the full 
range of s. H3 and H4 have been calculated by linear fit over the range 10 2 < s < 10 4 . One can 
observe that the quality of the fit improves in the central range of scales. 



H 


Hi 


Pi 


H 2 


P2 


H 3 


P3 


H 4 


Pa 


0.1 


0.1503 


0.9962 


0.1470 


0.9965 


0.1520 


0.9998 


0.1445 


0.9995 


0.2 


0.1973 


0.9984 


0.1939 


0.9984 


0.2005 


0.9999 


0.2028 


0.9999 


U.o 


0.2700 


0.9998 


0.2688 


0.9998 


U.Z 1 OO 


0.9999 


0.2809 




0.4 


0.3780 


0.9984 


0.3881 


0.9980 


0.3718 


0.9993 


0.3558 


0.9997 


0.5 


0.4467 


0.9994 


0.4513 


0.9994 


0.4514 


0.9997 


0.4656 


0.9993 


0.6 


0.5362 


0.9992 


0.5418 


0.9993 


0.5494 


0.9996 


0.5426 


0.9997 


0.7 


0.6129 


0.9996 


0.6115 


0.9997 


0.6349 


0.9998 


0.6877 


0.9992 


0.8 


0.7430 


0.9993 


0.7412 


0.9995 


0.7741 


0.9997 


0.7863 


0.9997 


0.9 


0.8645 


0.9990 


0.8747 


0.9991 


0.8650 


0.9999 


0.8775 


0.9999 


means 


of fractal descriptors. 


Generally, such a 


structure can 


properly reproduce 


complex 



systems whose heterogeneity is described by correlation decaying as a power law over space. 
This correlation corresponds to the intermediate behavior of real structure, as opposed to 
the fully uncorrelated exponential decay and the fully correlated constant decay exhibited 
respectively by the correlation function of ideal cases such as the perfect gas and the regular 
lattice. 
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